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ABSTRACT 

We describe a new numerical algorithm to obtain high-resolution simulated maps of the Cosmic Microwave 
Background (CMB), for a broad class of non-Gaussian models. The kind of non-Gaussianity we account for is 
based on the simple idea that the primordial gravitational potential is obtained by a non-linear but local mapping 
from an underlying Gaussian random field, as resulting from a variety of inflationary models. Our technique, 
which is based on a direct realization of the potential in spherical coordinates and fully accounts for the radiation 
transfer function, allows to simulate non-Gaussian CMB maps down to the Planck resolution (^ max ~ 3,000), with 
reasonable memory storage and computational time. 

Subject headings: Cosmic microwave background — Cosmology: theory — statistics — numerical simulations 



1. INTRODUCTION 

It has been only recently realized that a certain degree of non- 
Gaussianity is a general prediction of all inflationary models for 
the generation of primordial cosmological perturbations. Since 
the statistics of temperature anisotropies of the Cosmic Mi- 
crowave Background (CMB) directly probes that of the primor- 
dial fluctuations, any detection of non-Gaussianity in the CMB 
sky would provide a crucial test of the inflationary paradigm. 

The predicted degree of primordial non-Gaussianity is, how- 
ever, generally small; moreover, foreground contamination, in- 
strumental noise and many secondary sources of anisotropy also 
generate non-Gaussian features in the CMB. Cosmic variance 
further complicates the problem, as some non-Gaussianity un- 
avoidably arises from the uniqueness of the observed CMB sky 
(Scaramella & Vittorio 1991; Srednicki 1993; Gangui et al. 
1994; Gangui & Martin 2000; Komatsu & Spergel 2001). All 
these effects make the search for non-Gaussianity a challenging 
observational and statistical task. 

Because of these facts, the quest for primordial non-Gaussian 
signals in the CMB requires high-resolution all-sky measure- 
ments, such as those recently achieved by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) 3 satellite (e.g. Spergel et 
al. 2003) and those that the forthcoming Planck Surveyor 4 
will provide. At the same time one will need to develop op- 
timal statistical estimators specifically designed to search for 
non-Gaussianity of a given type. This makes extremely im- 
portant testing the power of different estimators on simulated 
all-sky CMB maps, with resolution comparable to that of the 
WMAP and Planck experiments. The aim of this paper is pre- 
cisely that of providing an efficient numerical algorithm able to 
produce simulated CMB maps for a broad class of physically 
motivated non-Gaussian models. Simulations of non-Gaussian 
maps of the CMB sky have been recently performed with a 
number of different techniques (Contaldi & Magueijo 2001; 
Vio et al. 2002; Martinez-Gonzalez et al. 2002; Komatsu et 
al. 2003). The kind of models that we consider (see also Ko- 
matsu et al. 2003) are based on the simple idea that the primor- 
dial (i.e. before being processed by the post-recombination ra- 
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diation transfer function) gravitational potential $(x) [actually 
Bardeen's gauge invariant variable $# (Bardeen 1980)] can be 
written as a local but non-linear mapping from a 'linear' ran- 
dom field $ L (x), which is Gaussian distributed with zero mean, 
namely 

$(x) = ^[$ L (x)]-(^[$ L ]} , (1) 

where the last term on the RHS, has been added to ensure that 
$ has zero mean, and only affects the mean temperature of the 
CMB. Models of this type were originally proposed by Coles 
and Barrow (1987), in connection with the statistics of CMB 
anisotropies in the large-angle (Sachs-Wolfe) limit. Moscardini 
et al. (1991) considered various models belonging to this class, 
as initial conditions for N-body simulations of large-scale struc- 
ture formation in a Cold Dark Matter (CDM) cosmology. Falk, 
Rangarajan and Srednicki (1993) showed that perturbation pro- 
duced during inflation, in single-field models, are characterized 
by a low non-Gaussianity level, which can be parametrized as 
(e.g. Wang & Kamionkowski 2000) 

$(x) = $ l (x) + /nl ($ l (x) 2 - <($ l ) 2 » , (2) 

where the dimensionless 'non-linearity' parameter /nl is much 
smaller than unity in their calculation, as a direct consequence 
of the tiny inflaton self-coupling. Similar conclusions were 
reached by Gangui et al. (1994), who used the techniques of 
stochastic inflation, to account for the back-reaction of field 
fluctuations on the background metric. Related analyses have 
been made by Salopek, Bond & Bardeen (1989), Salopek & 
Bond (1990, 1991) and Wang & Kamionkowski (2000), who 
accounted for possible features in the inflaton potential; Gupta 
et al. (2002), extended the stochastic calculation to 'warm- 
inflation' models, to evaluate the expected amplitude of non- 
Gaussianity. 

The key idea in this field is to adopt the model described 
by Eq.(2) as a prototype of the more general class of Eq.(l), 
once the functional T is Taylor expanded beyond the linear 
term, as one can always do, given that only a mild level of 
non-Gaussianity is allowed by the data. Verde et al. (2000) 
showed that the optimal strategy to constrain this kind of mod- 
els is to use (higher-order) statistics of the CMB, while Matar- 
rese, Verde & Jimenez (2000) and Verde et al. (2001) consid- 
ered this class of models in connection with the abundance of 
high-redshift objects, which can be mapped to rare events of the 
underlying PDF (see also Komatsu et al. 2003). 

Quite recently, there has been a burst of interest for non- 
Gaussian perturbations of this type. Different CMB datasets 
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have been analyzed, with a variety of statistical techniques (e.g. 
Komatsu et al. 2002; Cayon et al. 2003; Santos et al. 2003), 
with the aim of constraining the 'non-linearity' parameter / NL . 
The most stringent limit to date has been obtained by the WMAP 
team (Komatsu et al. 2003): -58 < / NL < 134 at 95 % cl. 
Komatsu & Spergel (2001) showed that the minimum value of 
|/nl| which can be in principle detected using the angular bis- 
pectrum, is around 20 for WMAP, 5 for Planck and 3 for an 
ideal experiment, owing to the intrinsic limitations caused by 
cosmic variance. Alternative strategies, based on the multivari- 
ate empirical distribution function of the spherical harmonics 
of a CMB map, have also been proposed (Hansen et al. 2002; 
Hansen, Marinucci & Vittorio 2003). 

On the theoretical side, a massive analytical effort has been 
recently made in order to obtain a quantitative prediction for the 
non-linearity parameter in slow-roll inflation models (Acqua- 
viva et al. 2002; Maldacena 2002). The resulting expression 
appears to be more complicated than previously thought, as the 
parameter / NL is replaced by a suitable convolution kernel /C, 
namely 

/ NL $ L (x) 2 ^ Jd 3 x 1 d 3 x 2 lC(x 1 -x,x 2 -x)<S> L (x l )<f> L (x 2 ). (3) 

Although the 'shape dependence' implied by the kernel (Ac- 
quaviva et al. 2002; Maldacena 2002) might be an interesting 
way to look for specific inflationary non-Gaussian signatures, 
the typical values assumed by JC are of order 10" 1 - 10" 2 , which 
makes the detection of non-Gaussianity produced in single-field 
slow-roll inflation a very challenging task! On the other hand, 
there are many physically motivated inflationary models, which 
can easily accommodate values of /nl as large as, or even much 
larger than unity. This is the case, for instance, of a large 
class of multi-field inflation models, which leads to either non- 
Gaussian isocurvature perturbations (Linde & Mukhanov 1997; 
Peebles 1997; Bucher & Zhu 1997) or cross-correlated and non- 
Gaussian adiabatic and isocurvature modes, as first realized by 
Bartolo, Matarrese & Riotto (2002) (see also Bernardeau & 
Uzan 2002). Other interesting possibilities include the so-called 
curvaton model, where the late time decay of a scalar field, 
belonging to the non-inflatonic sector of the theory, induces 
curvature perturbations (e.g. Mollerach 1990; Lyth & Wands 
2002) which are characterized by values of / N l > 1 (Lyth, 
Ungarelli & Wands 2003). The possible role of local fluctu- 
ations in the inflaton coupling to ordinary matter in driving 
non-Gaussianity has been recently analyzed (Dvali, Gruzinov 
& Zaldarriaga 2003; Kofman 2003; Zaldarriaga 2003). Non- 
Gaussianity caused by the effect of higher-dimension operators 
in the Lagrangian of single-field inflation models has also been 
considered (Creminelli 2003). 

The model considered here will be that of Eq.(2), which, be- 
sides the advantage of its simplicity, is very useful to constrain 
the non-Gaussianity level in terms of a single parameter / NL . 

The plan of the paper is as follows. In Section 2 we discuss 
the various steps which lead to the construction of our simula- 
tion algorithm. Section 3 contains our main results and some 
concluding remarks. 

2. OUTLINE OF THE METHOD 

The primordial potential fluctuations can be related to tem- 
perature anisotropies of the CMB by convolution with the radi- 
ation transfer function Ai(k), 
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FIG. 1 . — The coefficients A((r) evaluated for different values of £, assuming 
a ACDM cosmology (h = 0.65, Cl b = 0.05, Q c = 0.25, Q A = 0.7). We put r = cq, 
where r) is the conformal time. In the chosen model, crfy ~ \A.9Gpc and 
C7f* ~ 289 Mpc, where r/ and r/, are the present day and the recombination 
conformal time respectively. 



where the temperature fluctuation in the direction h has been 
expanded in spherical harmonics Yg m (n) with multipole coeffi- 
cients ai m , namely ^?(n) = ^2 lm ai m Yg m {n). The problem ana- 
lyzed here is how to simulate non-Gaussian temperature maps 
starting from the non-linear primordial gravitational potential 
described by Eq. (2). To this aim it is necessary to obtain the 
non-Gaussian part of the potential as a function of the Gaus- 
sian one, which can be generated by standard techniques, and 
finally account for the radiative transfer, in order to determine 
the Gaussian and non-Gaussian contributions to the tempera- 
ture harmonic coefficients ai m . At face value the simplest ap- 
proach seems to make realizations of the primordial potential in 
Fourier space and then replace it directly in equation (4). How- 
ever, this approach presents some technical drawbacks, as we 
are going to explain. In Fourier space, the non-Gaussian part 
of the potential is a convolution product of Gaussian fields (e.g. 
Komatsu & Spergel 2001): 

$ NL (k) = fa[£P- $ L (k+p)<I> L (p) , (5) 

J (2tt) j 

so, in principle we may obtain $ NL (k) by a numerical evalu- 
ation of the convolution, by starting from a Gaussian poten- 
tial field with given power-spectrum, generated on a Fourier- 
space grid. The computational problems in this procedure arise 
from the size of the simulation box. We have in fact to con- 
sider a box-size of the order of the present day cosmic horizon, 
Lbox = 2cr]o ~ 20 h~ x Gpc, (where 770 is the value of the confor- 
mal time today) and we need a resolution of about 20/T 1 Mpc 
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to accurately resolve the last scattering surface (e.g Komatsu et 
al. 2003; see also Figure 1): from this we find that, in order 
to generate the Gaussian part of the potential we need at least, 
1024 3 points on the grid. 

5 

So, the recipe for the generation of the non-Gaussian maps 
could be the following: i) generate a Gaussian potential in Fourier 
space, using 1024 3 grid-points; ii) go to real space by inverse 
Fast Fourier Transform (FFT), and square it to obtain the non- 
Gaussian contribution 6 ; Hi) transform back to Fourier space and 
convolve the harmonic transforms of the Gaussian and non- 
Gaussian parts [see Eq. (4)] with the radiation transfer function 
A e (k), to get the CMB multipoles. If we call <5>i m (k) the multi- 
poles of the harmonic expansion of the primordial gravitational 
potential, the CMB multipoles are in fact given by 
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where <$> £m (k) = ^Jk) + f NL ^(k) and, for each term, 
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The method we have just outlined can be easily improved 
by accounting for radiative transfer in real rather than Fourier 
space: to see this let us start from the harmonic transform of the 
primordial potential in real space and call it <f>£ m (r) [in analogy 
with the previous definitions we will also introduce the quanti- 
ties $\ m {r) and We can find the relation between any 
pair of &im(k) and &t m (r) (see the Appendix), namely 



2tt 2 



dkk 2 <$> lm (k)j t (kr) , 



and its inverse 



<f> em (k) = 4ir(i) e / drr 2 <f> em (r)j e (kr) 



(8) 



(9) 



where jg is the spherical Bessel function of order I. 

If we now replace Eq. (9) in Eq. (6), exchange the order of 
integrations and define 



Ai(r) =~ J dkk 2 Ai(k)je(kr) , 
we can finally write: 

a im = / drr 2 $ iln (r)A e (r) . 



(10) 



(11) 



The numerical advantages deriving from Eq. (11) are clear: 
if we generate CMB multipoles by this equation we can pre- 
compute and store the quantities Ag{r) for all the simulations 
characterized by the same cosmological model and we have to 
FFT only once, instead of going to real space and then back to 
Fourier space. The numerical calculation of the Ag{r) has been 
obtained by a modification of the CMBfast code (Seljak and 
Zaldarriaga 1996); the behaviour for some of these coefficients, 
evaluated for a ACDM cosmology, is shown in Figure 1 . 

To summarize: a numerically efficient technique for the gen- 
eration of non-Gaussian maps is to generate a Gaussian field on 
a grid in Fourier space, then FFT and square it to get the non- 
Gaussian part of the gravitational potential in real space and 
finally harmonic transform and convolve the two contributions 
with A^(r) to obtain Gaussian a\ m (r) and non-Gaussian ones, 
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FIG. 2. — The top panel shows the spherical Bessel function of order 150. 
The bottom panel shows the filters Wi(r, r[ ) as functions of r\, at fixed r, and 
three different values of I. We choose the spectral index n = 1 . Note that the 
filters become more and more peaked around r\ = r as ( increases. The Bessel 
functions oscillate very fast and have to be sampled in many more points than 
Wi to have an accurate numerical integration. For this reason the use of Wg 
instead of je drastically reduces the needed CPU time, as explained in the text. 



a^(r). For any value of /ml the CMB multipoles will then be 
given by ag m = a\ m +/ NL a^. 

Komatsu et al. (2003) (see also Komatsu, Spergel & Wandelt 
2003) used this method to generate 300 non-Gaussian maps at 
the WMAP resolution; their algorithm required about 3 hours 
on 1 processor of SGI Origin 300 and 7.5 GB of physical mem- 
ory to generate a single map. 

We developed an alternative algorithm which starts from the 
generation of the linear potential multipoles $^ m (r), while we 
use the same technique as Komatsu et al. (2003) to account for 
the radiative transfer. We believe that the main advantages of 
the method we are going to describe w.r.t. the one just outlined 
can be the following: i) much less memory requirements, ii) 
a decrease of CPU time, Hi) the possibility to more accurately 
resolve the last-scattering surface. 

As we are going to see, we will obtain the multipole coef- 
ficients &\ m {r) directly, i.e. without passing through the gen- 
eration of the potential on a cubic grid. So we won't need 
to convert from Cartesian to spherical coordinates, in order to 
harmonic transform $ L (r) and, even more important, we won't 

5 Note that at any stage one can replace the square operation by whatever local 
functional of with no extra computational complications. 

6 Note that one can replace the square operation by whatever local functional 
of <3>l, with no extra computational complications. 
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need to store big arrays containing the values of $ L (r), calcu- 
lated on the grid points. As the <&\ m (r) are Gaussian variables, 
to generate them we only need to know their correlation func- 
tion. We have (see the Appendix): 

(12) 

where Pg,(k) is the primordial (i.e. unprocesssed by the radia- 
tion transfer function) power-spectrum of the gravitational po- 
tential, and S'j is Kronecker's delta. Since the multipoles <& L m 
are correlated in real space, they cannot be generated directly. 
One possibility is to generate $\ m (k) first, as they are indepen- 
dent Gaussian variables (see the Appendix, for an explicit cal- 
culation), 

<*/ lBl (*i)*&*(*i)> = 8^ 3 ^V(£i > (") 

(5 D being the Dirac delta function) and then transform <&i m (k) 
into <$>tm(r), by Eq. (8). 

As in this last equation, I and m are fixed and only k varies 
for each integration, it is clear that this method does not require 
large memory: in fact, we never need to store all the values of 
&e m (k) and &g m (r) at the same time in physical memory. 
Unfortunately, the spherical Bessel transform which allows the 
passage from &e m (k) to &i m (r) requires a lot of CPU time, as the 
jt(kr) are rapidly oscillating functions which have to be sam- 
pled in many points to ensure accurate numerical integration. 
The time needed to simulate one sky-map with this method is 
around 20 hours on a Compaq server DS20E (500 Mhz), for 
Imwc = 500, and scales roughly 7 as I 2 . However, it is clear from 
the previous argument that the code could be greatly speeded up 
by avoiding the numerical integration of Eq. (8). This result can 
be achieved as follows: start from independent complex Gaus- 
sian variables nt m (r), characterized by the correlation function 

(n tini (n)n* i2m2 (r 2 )) = ^"^ C > < 14 ) 

as shown in the Appendix, it is now possible to recover our 
Gaussian variables <&\ m (r), with the right correlation properties, 
by convolution, 

= J dn rlnUnW^n) , (15) 

where the filter functions Wg(r, r\) are defined as 

W e (r,n) = | J dkk 2 y/PrffiUmMkr!) . (16) 

In this work we will only consider models with a scale-free 
primordial density power-spectrum of the form P(k)=Ak", and 
we take n = 1 , in agreement with observational data (e.g. Spergel 
et al. 2003). In this case, when r is fixed, the functions We(r, r\) 
are smooth and differ from zero only in a narrow region around 
r = r x (see Figure 2). So, in order to accurately sample for 
their integration in Eq. (15) we only need to compute them in 
a few points, which considerably reduces the CPU time needed 
for the generation of ^\ m (r). The CPU time needed for the gen- 
eration of a single sky-map at the WMAP resolution now goes 
down to less than one hour. In other words, by the latter method 
we work directly in real space, starting from white-noise coef- 
ficients and convolving them with suitable filters to produce the 

7 This is because the total number of multipoles, and consequently the number 
of transforms to be calculated, is ~ I 2 



right correlation properties of the $\ m {r). This technique for 
the generation of correlated random variables is well known in 
the literature (e.g. Salmon 1996; see also Peebles 1983). It is 
particularly efficient in our case, because, as it appears from 
Figure 2 the We(r, r{) get more and more peaked around r= r\ 
as £ increases. 

3. RESULTS AND CONCLUSIONS 

Let us summarize the various steps which define our numer- 
ical algorithm. 

• Simulate the white-noise Gaussian fields ne m (r) directly 
in real space. 

• Convolve them with the (pre-computed) filter functions 
Wg(r, ri), to obtain the linear potential multipoles $\ m (r). 

• Compute the linear potential $ L (r) = Y,i m ®e m (r)Ye m (r) 
and square it to obtain (modulo / NL ) the non-Gaussian 
contribution to the gravitational potential, directly in spher- 
ical coordinates. 

• Harmonic transform 8 <£> NL (r) to get &^(r). 

• Convolve <$> lm {r) = &i m {r) + /nl^W with the (pre- 
computed) real-space radiation transfer functions Ae(r). 

The form of Ai(r) suggests to generate the potential coeffi- 
cients $£ m (r) for non-uniformly spaced values of the radial co- 
ordinate. In the final convolution, the real-space radiation trans- 
fer functions are in fact different from zero only near last scat- 
tering (0. lr/* < rj < 277*, see Figure 1 ; 77* is the conformal time, 
correspondingly = cry*), with a further contribution, coming 
from low redshifts and multipoles (£ < 50, r < 5000 /T 1 Mpc, 
see the lower panel of Figure 1), due to the late-ISW effect. 
Therefore, we generated the multipoles &( m (r) only for I and r 
belonging to the previously defined intervals 9 . In this way we 
can increase the radial resolution in the selected regions and, at 
the same time, reduce the total number of points in our spheri- 
cal grid (i.e. reduce the total number of evaluated $i m (r)), with 
a consistent reduction of the total CPU time needed. More- 
over, since at low redshifts we need to consider only low mul- 
tipoles, we may in this case reduce the pixel resolution when 
using HEALPix to perform the harmonic transforms, with fur- 
ther reduction of CPU time (this is because the value of the 
parameter nside in HEALPix, which is directly related to the 
pixel resolution, has to be chosen in such a way as to verify 
the relation l max < 3 x nside, and the calculation time for the 
harmonic transforms scales as (nside) 3 ). 

In other words, our multipole approach enables to adapt the 
simulation grid to the optimal sampling of the physically rele- 
vant radial intervals (i.e. points where the radiative transfer is 

8 Harmonic transforms have been performed using the HEALPix pack- 
age [http://www.eso.org/science/healpix/ see also (Gorski, Hivon & Wandelt 
1999)]. 

9 Actually, at low redshifts we generate ^ m (r) for I < 100, even if we are 
ultimately interested only in the multipoles with i < 50; this is necessary as 
high-order Gaussian multipoles contribute to lower-order non-Gaussian ones in 
the harmonic transforms. 



5 



3x10* 
2x10* 
'0' 


3x10* 
2x10* 
•0' 




' : f NL =500 - 


- f NL =1000 - 


' : f NL =3ooo : 


'- ... f NL -5ooo : 

f. 



-400 -200 200 400 -400 -200 200 400 

AT (/xK) AT (fiK) 

FIG. 3. — One-point PDF of temperature fluctuations obtained from a single 
simulated map with l mar = 3000. Different panels show a comparison between 
the Gaussian realization (/nl = 0) and the corresponding non-Gaussian one, 
for different values of /nl- Beam smearing has not been included here. 



non-zero) and to reduce the CPU time needed for the genera- 
tion of (r). Moreover, we found it useful to split the high 
redshift contributions (SW effect, acoustic oscillations) from 
the low redshift contributions (late-ISW effect): in this way the 
simulation of the latter effects can be greatly speeded up, as 
they affect only the lower multipoles of the final map. 

We have applied our algorithm to obtain some simulated tem- 
perature maps with different values of £ max , starting from a 
ACDM model with primordial spectral index n = 1 . The radial 
coordinate has been discretized with a step of 7 .5 h~ l Mpc in the 
regions where Ai(r) =/0. In this way we generated 100 poten- 
tial multipoles in the radial interval (rn— 2r*) < r < (ro-O.lr*), 
for each £,m, 2 < £ < / max , < m < £, and we generated about 
600 potential multipoles with r < 5000 h Mpc, for each 2 < 
£ < 100, 0<m<£. 

To run our simulations we need 200 Mb of physical memory 
when ^ max = 3000; the final CPU time is reported in Table (1), 
where we also show the time needed for the harmonic analysis, 
at large r only, as this is the most time-consuming part of our 
algorithm. 

In Figure 3 we plot the PDF of temperature fluctuations, ob- 
tained from an original simulation with £ max = 3000, mainly 
for a comparison with the results contained in (Komatsu et al. 
2003). In Figures 4,5 we show maps obtained from the same 
simulation, smoothed with different Gaussian beams at the res- 
olution of COBE DMR, WMAP and Planck. 



The main purpose of obtaining simulated non-Gaussian CMB 
maps is that of providing a test-bed for the power of different 
statistical estimators, specifically designed to search for non- 
Gaussianity, in detecting the kind and the level of non-Gaussianity 
which is present in the maps. To this aim, all the possible ex- 
tra effects, such as foreground contamination and instrumental 
noise of a specific experiment, which appear in the statistical 
analysis of real datasets should be properly taken into account. 
These aspects will be the analyzed in a forthcoming paper. 

As already mentioned, there are various non-primordial cos- 
mological effects which also imply a deviation from the Gaus- 
sian behaviour. These include systematic effects as well as 
some secondary sources of anisotropy. Among the future de- 
velopments of the present algorithm is the possibility to include 
suitable kernels to account for some of these effects in the sim- 
ulated CMB maps. This issue will be discussed elsewhere. 

APPENDIX 
DERIVATION OF USEFUL FORMULAE 
Relation between &£ m (k) and §i m (r) 
We start from the Fourier transform 

$ W= / 7TT3^ k ) e "' k ' r (A V 
J (2ttY 

and we Rayleigh expand exp(-/k ■ r) in spherical harmonics, 
expHk-r] = 4tt £ ^(/)- £ T, m (k)F; m (f )j t (kr) ; (A2) 

£ m 

we find: 

m = EEH)'*7 m (*) / dkk 2 j e (kr) x 

x y<i%y^(£)$(k). (A3) 

Recalling that: 

*im(r) = Jdn f ^(r)Y; m (r) (A4) 

we can replace the previous expression for $(r) in the last for- 
mula [definition of Q( m (r)] and expand also ^(k) in spherical 
harmonics. We obtain 

®Ur) = ^2 E E("'') £l / dil f Y; m (r)Yl mi m x 

e,e 2 »«i'«2 

x Jdkk 2 j ei (kr)<^ e2m2 (k) yrfO £ ^ imi (k)y 4m2 (£)(A5) 

It is now easy to see, by the orthonormality of spherical har- 
monics, that this reduces to Eq. (8). 



^max 


CPU time 


Harmonic transform 


nside 


300 


8 min 


6 min 


128 


500 
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32 min 


256 


750 


1 h 8 min 
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256 


1500 


8 h 47 min 


8 h 20 min 


512 


3000 


68 h 45 min 


67 h 15 min 


1024 



Table 1 

CPU time needed to simulate a single map 
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FIG. 4. — The left-hand side of this panel shows the same Gaussian realization, smoothed by three different beams. From top to bottom, the FWHM of the beams 
is 1° (COBE), 13 arcmin (WMAP) and 5 arcmin (Planck). The right-hand side shows three corresponding non-Gaussian realizations, obtained from the Gaussian 
one by adding a non-linear coupling parameter /nl = 3000 (such a high value of /nl is chosen to make the non-Gaussian effects visible by eye). The model is a 
ACDM with primordial spectral index n=l. The primordial fluctuations have been COfiE-normalized by CMBfast. 
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<mapG_Planck.jpg> 
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FIG. 5. — In this panel we show the same Gaussian maps of Fig. 4 on the left-hand side, while on the right-hand side we take /jvl = -3000. The remaining 
parameters are the same as in Fig. 4 



Correlation function of§\ m (k) 
Using the definition of $\ m {k), we can write 



<*fc Bl (*i)*Wfe)> = J dSkJSlt, (<J> L (k!)<i> L *(k 2 )) x 

xy^kDF^fe). (A6) 

From this we obtain: 

($Wfci)$£ m2 (fc 2 )) = ( 2 3p *(*i) / «ttVKV D (ki-k2) x 

xy,* mi to 2m2 (k 2 ). (A7) 

Where we have used: 

($ L (k 1 )$^ (k2 )) = (2 7 r) 3 P tI) (fe 1 )^(k 1 -k 2 ) . (A8) 
Recalling now the integral representation of Dirac's delta func- 
tion, 

Akx - k 2 ) = — fd 3 r exp[i(ki - k 2 ) -r] (A9) 
we can Rayleigh expand to obtain: 

l\l- L t 1 ,l i m\m2m- i mt, 

x Jdn- r Y ijmj (r)Y; im4 (r) J dS^Y^foY^fa) x 

x y jf> fe2 y, 2 „ 12 (k 2 )y, 4m4 (k 2 ) . (aio) 

By the orthonormality of spherical harmonics this reduces to 

(^; mi (A:i)$£ m2 (fc 2 )) = 167r 2 P$(A:i)^ lfe (5 mim2 x 

xjdrr 2 j ei (kir)j i2 (k 2 r), (All) 

from which, recalling that 

/■OO 1 

jf drr 2 j e (k l r)j e (k 2 r) = --^8 D (k l -k 2 ), (A12) 
we finally obtain 

{$i mi (kx)&£ m2 (k 2 )) = ^^-s^-^spz: . (AO) 

Correlation function of <&\ m {r) 
According to the Wiener- Khintchine theorem, 

J (2tt) 3 

where r = |ri-r 2 |, and P$(k) are the correlation function 
and the power spectrum of primordial potential fluctuations, re- 
spectively. Expanding once again the exponential in spherical 
harmonics, we obtain 

£* (r) = - Y im(h)Y: m (h) / dkk 2 P^k)Mkn)Mkr 2 ) . 

77 1=0 m=-i J 

(A15) 

By the definition of $t m (r), we can write 
(^: mi (n)^i m2 (r 2 )) = J dn h dn f j; imi (h)Y l2m (r 2 )te(r) . 

(A16) 

The final result is obtained by substituting by its expres- 
sion (A15) and using the orthonormality of spherical harmon- 
ics, 

(*£ mi (r!)#Uta)> = \^\ 8 Z J dk&PtWfaQrOMkri) ■ 

(A17) 



Derivation ofWp(r,r\) 

Consider a Gaussian random field $ L (r) with correlation func- 
tion £$(r). If P<s>(k) is the corresponding power spectrum, then 
$(r) can be obtained as: 

$ L (r) = J d\Mn) J ^VP^ky k (c ~ Cl) , (A18) 

where we have introduced the white-noise field n(r{): 

(«(r 1 )«(r 2 )) = ( 5 D (r 1 -r 2 ). (A19) 
To find a similar expression for $\ m {r) we define the quantities 

W(r-n) = / ^y/IwW*-™ ; (A20) 

J (27T) 3 

now, as usual, we expand both $(ri) and W(r-ri) in spheri- 
cal harmonics and account for the orthonormality relations, to 
finally obtain 



<&\jr) = I dri r\n lm {n)W t {r,n) 



(A21) 



W e (r, n) = ^J dkk 2 y/P*(k)jt(kr)je(kn) (A22) 



and 



(n^(n)nl m2 (r 2 ))= 5D{r l 2 ^ 8%6% 



(A23) 
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D. Marinucci for useful discussions. 
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